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ABSTRACT 

We present numerical results on the properties of young binary and multiple stel- 
lar systems. Our analysis is based on a series of SPH + A^-body simulations of the 
fragmentation of small molecular clouds, that fully resolve the opacity limit for frag- 
mentation. These simulations demonstrate that multiple star formation is a major 
channel for star formation in turbulent flows. We have produced a statistically signif- 
icant number of stable multiple systems, with components separations in the range 
~ 1 — 10"^ AU. At the end of the hydrodynamic stage (0.5 Myr) we find that w 60% of 
stars and brown dwarfs are members of multiples systems, with about a third of these 
being low mass, weakly bound outliers in wide eccentric orbits. Our results imply that 
in the stellar regime most stars are in multiples (« 80%) and that this fraction is 
an increasing function of primary mass. After iV-body integration to 10.5 Myr, the 
percentage of bound objects has dropped to about 40%, this decrease arising mostly 
from very low mass stars and brown dwarfs that have been released into the field. 
Brown dwarfs are never found to be very close companions to stars (the brown dwarf 
desert at very small separations) ^ but one case exists of a brown dwarf companion at 
intermediate separations (10 AU). Our simulations can accommodate the existence 
of brown dwarf companions at large separations, but only if the primaries of these 
systems are themselves multiples. 

We have compared the outcome of our simulations with the properties of real 
stellar systems as deduced from the IR colour-magnitude diagram of the Praesepe 
cluster and from spectroscopic and high-resolution imaging surveys of young clusters 
and the field. We find that the spread of the observed main sequence of Praesepe in 
the 0.4 — IM0 range appears to require that stars are indeed commonly assembled into 
high-order multiple systems. Similarly, observational results from Taurus and p Ophi- 
uchus, or moving groups such a TW Hydrae and MBM 12, suggest that companion 
frequencies in young systems can be indeed as high as we predict. The comparison with 
observational data also illustrates two problems with the simulation results. Firstly, 
low mass ratio (g < 0.2) binaries are not produced by our models, in conflict with 
both the Praesepe colour magnitude diagram and independent evidence from field 
binary surveys. Secondly, very low mass stars and brown dwarf binaries appear to be 
considerably under-produced by our simulations. 
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1 INTRODUCTION 

The formation of binary and multiple systems provides one 
of the most exacting areas in which star formation the- 
ory can be compared with observational data. Stars are 
known to have a binary frequency in excess of 50%, both 
in the field (Duquennoy & Mayor 1991, henceforth DM91; 
Fisher & Marcy 1992; Halbwachs et al. 2003) and in clus- 
ters (Kohler & Leinert 1998; Patience et al. 1998; Math- 
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ieu et al. 2000; Bouvier et al. 2001). For pre-main-sequence 
stars this frequency seems to be even higher (Reipurth & 
Zinnecker 1993; Simon et al. 1995; Duchene 1999; Reipurth 
2000). Thus, understanding the formation of multiple stars 
becomes necessary if we are to understand star formation 
in general. Recently, Bate, Bonnell & Bromm (2003; hence- 
forth BBB) have shown that multiple stars are a natural 
byproduct of the collapse and fragmentation of turbulent 
molecular clouds. In particular, close binary stars can be 
indirectly formed out of the fragmentation of a molecular 
cloud core by means of dynamical interactions with other 
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stellar objects in the vicinity, as well as interactions with cir- 
cumbinary discs (Bate, Bonnell & Bromm 2002b). However, 
this sort of calculation (that fully resolve fragmentation in 
a cloud of e.g. 50 Mq of gas in BBB's case) is very demand- 
ing computationally, and a wide range of initial conditions 
cannot easily be explored. In addition, such calculations re- 
main inconclusive about the properties of multiple stars (in 
particular the existence of long lived multiples and the pro- 
duction of wide binaries) since the evolution of the resulting 
stellar systems is not followed until decay to a stable hier- 
archical configuration is achieved. For example, out of the 
50 objects formed in the BBB simulation, 18 are contained 
in non-hierarchical groupings containing 11 and 7 members. 
Since, in this simulation, only the more compact groupings 
have had time to decay to a stable configuration, BBB un- 
derstandably focused their attention on the implications of 
their work for close binary formation. 

In this paper we have taken an alternative approach 
which is a natural complement to the BBB simulations. 
Through the modelling of less massive clouds over longer 
timescales, we have been able to follow wide binary for- 
mation and the production of long lived multiple systems. 
We find that by modelling an ensemble of isolated cores of 
mass 5 Mq, we can improve the number of stars formed per 
CPU hour by a factor 7 compared with the BBB calcula- 
tion, which comprised a total mass of 50 M©. This economy 
stems from the fact that by focusing on individual dense 
cores, we dispense with the computational expense of fol- 
lowing the diffuse gas in the BBB simulation. We remark, 
however, as a caveat, that we find that the formation of mul- 
tiples proceeds hierarchically, with structures on any scale 
being progressively modified by interactions on larger scales. 
Thus in these models with scale free turbulence, any upper 
mass cut-off (i.e. the finite mass of the system modelled) 
may have some influence on the multiple systems produced. 
By comparing our results with those of BBB, we are able to 
assess how sensitive our results are to the total mass of the 
region simulated. 

The computational economies that we gain through the 
simulation of less massive cores, allows us to instead con- 
centrate our efforts on the longer timescale integration of 
the multiple systems produced. We wish to understand how 
our results are aifected if we throw away the gas component 
at a certain evolutionary stage (when 60% of the initial gas 
has been accreted) and thereafter evolve the system as an 
A'^-body ensemble. We follow the stellar dynamical A'^-body 
problem for 10 Myr (this typically corresponding to 10* or- 
bits of multiple systems at the median separation) , using the 
nbodyI code, by Aarseth (Aarscth 1963). Another advan- 
tage that we reap from this less computationally expensive 
approach is that we gain improved statistics (i.e. we form 
« 150 stars and brown dwarfs, compared with the « 50 of 
BBB), and from this dataset are able to extract reasonable 
statistics, in a manner similar to Sterzik & Durisen (1998; 
2003) for the purely A^-body case. We stress that our simu- 
lations share the property of the BBB simulations that they 
resolve the opacity limit for fragmentation (Low & Lynden- 
Bell 1976; Rees 1976) and that, assuming that fragmentation 
does not occur at densities greater than those at which the 
gas becomes opaque to infrared radiation, these calculations 
are able to model the formation of all the stars and brown 
dwarfs that, under the initial conditions imposed, can be 



produced. Our spatial resolution limit for binaries allows us 
to study a wide range of separations, and the particle num- 
bers we employ allow us to model accretion discs around 
the protostars which are as long lived as those modelled by 
BBB. 

Our study has four major findings: 

• The incidence of multiple systems (with N > 2) is high 
(8 out of 18 N > 1 systems aXt = 10.5 Myr) and we therefore 
test whether our results are compatible with the photometric 
width of the main sequence in young clusters. The compan- 
ion frequency varies significantly during the pure AT-body 
evolution of the systems. 

• Although we produce a variety of triples, quadruples 
and higher order multiples, they tend to follow a charac- 
teristic pattern of internal mass distribution. This involves 
each binary pair having a mass ratio that docs not deviate 
strongly from unity, and likewise, for quadruples, the mass 
ratio of the total mass in each binary is not greatly different 
from unity (i.e. mass ratios in the range 0.5 — 1 are favoured 
in each case) . The exception to this is that multiple systems 
are commonly orbited by low mass outliers on eccentric or- 
bits, which are the result of incomplete ejections of low mass 
objects from multiple systems. We discuss how deep images 
of the environs of multiple systems, and spectroscopic stud- 
ies of the primaries of extreme mass ratio binaries, may be 
used to test this prediction. 

• The multiplicity fraction is an increasing function of 
primary mass. The brown dwarf desert at very small sepa- 
rations is reproduced by our models. 

• We confirm the result of previous simulations that 
brown dwarf binaries indeed appear to be under-produced 
by turbulent fragmentation calculations. 

The structure of this paper is as follows. In Section 2 the 
computational method and initial conditions applied to our 

models arc described. The results on multiple stars are given 
in Section 3. In Section 4 we perform a detailed comparison 
of our results with available observational data, and suggest 
future experiments. Our conclusions are given in Section 5. 



2 COMPUTATIONAL METHOD 

The results presented in this paper were obtained from the 
same simulations introduced in Delgado-Donate, Clarke & 
Bate (2004; henceforth DCB04). A detailed description of 
the SPH code that was used to perform those calculations 
as well as the initial conditions imposed is given there. Here 
we briefly summarise the computational method. 

The calculations were performed using a 3D hybrid SPH 
A'-body code, with variable smoothing lengths. The SPH 
equations are solved using a second order Runge-Kutta- 
Fehlberg integrator with individual timesteps for each par- 
ticle (Bate, Bonnell & Price 1995). We use the standard 
form of artificial viscosity (Monaghan & Gingold 1983) with 
strength parameters = 1 and f3v = 2. 

2.1 Equation of state 

We have assumed that the gas becomes optically thick when 
the density p reaches a critical value pc = 10~^^g cm~^. This 
density defines the so-called opacity limit for fragmentation. 
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which sets a minimum fragment mass of a few Jupiter masses 
(Low & Lynden-Bell 1976; Boss 1988; EBB). 

To model the opacity limit for fragmentation without 
performing full radiative transfer, we use an equation of 
state given by p = K , where p is the pressure and K is 
a measure of the entropy of the gas. The value of rj changes 
with density as: 



n = 



{ 1, P < 10-^^ g cm-3, 
\ 7/5, p > 10"^^ g cm-^ 



(1) 



The gas is assumed to consist of pure molecular hydrogen {fj, 
= 2), and the value of K is such that when the gas is isother- 
mal, K = Cs, with the sound speed Cs = 1.85 x 10 cms^ at 
T = 10 K. The pressure is continuous when the value of r) 
changes. 

2.2 Sink particles 



A sink particle is inserted when the central density of a 
fragment exceeds pa = lO^^^g cm^'^, well above the critical 
density pc- Sink particles are point masses with an accretion 
radius, so that any gas particle that falls into it and is bound 
to the point mass is accreted. In the present calculations, 
the accretion radius -Rsink is constant and equal to 5 AU. 
Therefore, discs around sink particles will be resolved only 
if their radii <; 10 AU. Sink particles interact with the gas 
only via gravity and accretion. 

The gravitational acceleration between two sink parti- 
cles is Newtonian for r > 4 AU. but is smoothed within 
this radius using spline softening (Benz 1990). The maxi- 
mum acceleration occurs at r ~ 1 AU; therefore, this is the 
minimum binary separation that can be resolved. 



2.3 Initial Conditions 

We have performed 10 different calculations of the fragmen- 
tation of a small-scale, turbulent molecular cloud, each of 
them under almost exactly the same initial conditions. Each 
cloud core is spherical, has a mass of 5 Mq , a radius of « 10* 
AU, and an initial uniform density of 10~^*g cm^"^. At the 
initial temperature of 10 K, the mean thermal Jeans mass is 
0.5 M0, i.e. the Jeans number of the cloud is 10. The global 
free-fall time of the cloud te is « 10* yr. 

We have imposed an initial supersonic turbulent veloc- 
ity field on the gas, in the same manner as Ostriker, Stone 
& Gammie (2001) and BBB. We generate a divergence- 
free random Gaussian velocity field with a power spectrum 
P{k) oc , where k is the wavenumber and a is the power 
index, which we have set to —3 in half of the simulations 
and —5 in the other half. The velocity field is normalised 
so that initially it is in equipartition with the gravitational 
potential energy of the cloud core. 

These simulations do not include magnetic fields, as we 
have tried to isolate a particular hydrodynamical fragmen- 
tation problem to characterise the properties of the resulting 
stellar systems. We have not included in our models any me- 
chanical or radiative feedback mechanism. This may be an 
appropriate choice, since the maximum stellar mass in these 
simulations does not exceed 1 M0 , whereas the most power- 
ful winds and photoionisation fronts in star-forming regions 
are produced by much more massive stars. 



2.4 Resolution 

The local Jeans mass must be resolved throughout the cal- 
culation (Bate & Burkert 1997; Truelove ct al. 1997; Whit- 
worth 1998), otherwise some of the fragmentation might be 
artificially enhanced or suppressed. In order to model a 5 M© 
cloud core with critical density pc we need to use 3.5 x 10® 
particles. 

Each of the hydrodynamic calculations that are dis- 
cussed in this paper required « 4000 CPU hours on the 
SGI Origin 3800 Computer of the United Kingdom Astro- 
physical Fluids Facility (UKAFF). 



2.5 AT-body calculations 

The hydrodynamic evolution of each cloud is followed until 
60% of the initial gas particles are accreted. At this point, 
the remaining gas is removed and thereafter the system is 
evolved as an A'^-body ensemble for 10 Myr, using the numer- 
ical code NBODYl, by Aarseth (Aarseth 1963; see Aarseth 
1999 for an updated description of the code). nbodyI is a 
simple A''-body algorithm which computes the gravitational 
forces between particles by direct summation. The numeri- 
cal scheme is based upon the expansion of the gravitational 
force in a fourth-order polynomial with divided differences. 
Individual timesteps are also implemented. The Newtonian 
potential is softened in a manner that mimics a Plummer 
sphere. 

NBODYl lacks any regularization scheme for the treat- 
ment of close binaries, in contrast with subsequent more 
sophisticated nbody codes. Close binaries can nevertheless 
be accurately integrated provided that a sufficiently small 
softening length and 77 parameter (which controls the size of 
the timesteps) is chosen. We took softening lengths of 0.1 x 
the softening length of the SPH simulations, and r; parame- 
ters smaller than 10~^, and found that energy and angular 
momentum were conserved to an accuracy better than one 
part in one million. Each AT-body simulation required only 
« 50 CPU hours to be completed. 



3 THE PROPERTIES OF THE MULTIPLE 
STARS 

The hydrodynamical evolution of the cloud produces shocks 
which decrease the turbulent kinetic energy that initially 
supported the cloud. In parts of the system, gravity begins 
to dominate and dense self-gravitating cores form and col- 
lapse. These dense cores are the sites where the formation 
of stars and brown dwarfs occurs. The turbulence decays 
on the dynamical timescale of the cloud (as found by Mac 
Low et al. 1998; Stone, Ostriker & Gammie 1998; and BBB, 
among others), and star formation begins just after 1 to 1.5 
global free-fall times ta- As mentioned before, the hydrody- 
namical calculations were stopped when 60% of the gas has 
been accreted. In terms of ts this means that, for the a = — 3 
and a = — 5 calculations (henceforth a3 and a5), we follow 
the evolution of the cloud for « Atn and 5. Sift , respectively 
(i.e. an average of « 0.5 Myr). Altogether, the calculations 
produce 145 stars and brown dwarfs. An analysis of the de- 
pendence of the statistical properties of the resulting stars 
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Table 1. Multiple systems configurations. Each main row shows the results obtained for a given simulation: the first 5 correspond to 
the a = — 3 case, and the following 5 to the a = — 3 case. The last row shows the results for all the simulations combined. For each 
simulation, every column except the first two are divided in two rows: the upper one refers to the results at t = 0.5 Myr, i.e., at the end 
of the hydrodynamic calculation; the bottom row correspond to the results at t = 10.5 Myr, i.e., at the end of the A^-body integration. 

The Configuration column provides a symbolic representation of the hierarchical distribution of the stars and brown dwarfs making up 
a given multiple. The numbers can be interpreted as arbitrary flags assigned to each of the components of a given multiple (in fact, 
object 1 is simply the first star to have formed in that calculation, object 2 the second and so forth). Two numbers linked by a dash 
and between parentheses represent a binary system. Subsequent dashes denote further bonds with other objects; and square and curly 
brackets enclose those objects that belong to sub-systeuis with higher order in the hierarchy. For example, in simulation a?tA there is 
one multiple system which consists of a binary (3-5), orbiting a single (dubbed 0), thus making the triple [0-(3-5)]. This triple system is 
bound to another binary (1-2), thus forming a quintuple system {[0-(3-5)]-(l-2)}. Additional objects enclosed in parentheses represent 
outliers: individual very low mass stars or brown dwarfs also bound to the multiple but moving in eccentric orbits at large separation 
from the CoM of the system. A plus sign is used to highlight that a given simulation produced two or three mutually unbound multiples. 



Simulation N^^^ 



Configuration'^ 



ScpSut [AU] 



q;3A/ 



16 



7 (((({[0 - (3 - 5)] - (1 - 2)} - 14) - 15) - 4) - 6) 
10 ({[0- (3 -5)] -(1-2)} -15) 



10* 

4 X 10^ 



a3B 



21 



13 ((({[(1 - 3) - 7] - (0 - 2)} - 18) - 6) - 10) 
16 {[(1 - 3) - 7] - (0 - 2)} 



103 

10^ 



a3C 



17 



11 
(3,2,2) 



6 (((([{[(0 - 2) - 1] - (3 - 4)} - (6 - 7)] - 10) - 13) - 12) - 5) 
10 [(0-2)-l] + (3-4) -I- (6-7) 



10^ 
102 



a3D 



14 



10 
(3,2) 



4 (((({[(0 - 1) - (11 - 10)] - (3 - 2)} - 6) - 12) - 8) - 13) 

9 [(0 - 1) - 6] + (3 - 10) 



5 X 10-^ 
10 



a3E 



a5A 



17 



8 ((({[(1 - 2) - (7 - 13)] - (0 - 1)} - 16) - 5) - 11) 
11 {[1 - 2) - (V- i:-!)] - (0 - 1)} 



12 



(5,2) 
(4,2) 



6 



([(0-1) -(7 -2)] -4) + (8-9) 
[(0 - 1) - (7 - 2)] + (8-9) 



7 X 10" 

7 X 



2 X 10^ 
10^ 



a5B 



6 

(2,2) 



1 (([(1-2) -(0-4)] -5) -6) 
3 (1-2) -I- (0-4) 



2 X lO'^ 
10 



a5C 



11 



6 ({[(0-5) -2] -3} -9) 

7 ([(0 - 5) - 2] - 9) 



10^ 

3 X 10* 



a5D 



16 



9 

(2,2) 



7 ((({[(4 - 0) - 5] - [(1 - 9) - 14]} - 15) - 11) 
12 (1-9) + (4 - 0) 



10* 
10 



a5E 



All 



14 (6,2,4) 
(4,2,2) 



2 ({[(5 - 6) - (4 - 8)] - 12} - 9) + (13-10) 
6 [(1 - 2) - (0 - 3)] + (5-6) + (13-10) 



[(1 - 2) - (0 - 3)] 



86 : 59 
5.') : 



40% , 20% ] + 40% 
:«i%. . 2%. - ()2%. 



30 



4 X 10^ 
2 X 10^ 



" Ntot : total number of stars and brown dwarfs formed in a given calculation. 

Nm; number of stars/brown dwarfs bound to a given multiple system; Nsi number of singles. 

In the last row Configuration gives: % of inner companions , % of outer companions (as defined in below) -|- % of unbound objects 
Nout: number of bound objects with high eccentricities and separations larger than 1000 AU from the innermost component of the 
multiple system. 

" Sepout: Separation between innermost and outermost components of the multiple, in AU. 
f The [0-(3-5)] triple remains unstable at 10.5 Myr. 



and brown dwarfs (e.g. the IMF) on the different initial con- 
ditions imposed is presented in DCB04. Here, we focus on 

the formation of multiple stellar systems and their internal 
structure (internal distribution of masses and separations). 

3.1 Internal structure 

Both sets of initial conditions result in the formation of a 
large number of multiple stars. The statistical properties of 
the multiple systems are not very sensitive to the slope of the 
initial turbulent power spectrum (DCB04), implying that 



they are more sensitive to the dynamical and competitive 
accretion processes within each mini-cluster rather than the 
large scale morphology of the cloud. Henceforth we analyse 
the combined dataset for the a3 and ab runs. 

Independently of whether the gas is initially dominated 
by large- (a5) or small- (a3) scale turbulent motions, a 
series of localised, dense pockets of gas are formed within 
the cloud. A pressure-supported object forms within each of 
these dense clumps, first in isolation but soon surrounded by 
an accretion disc. Initially the mass of the disc is compaxa^ 
ble to, and often greater than the mass of the central object. 
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Thus, the disc is prone to the appearance of gravitational 
instabihties which, in most cases, result in the fragmenta- 
tion of the disc into one or more protostellar objects (Bon- 
nell 1994; Bonnell & Bate 1994; Burkert, Bate & Boden- 
hcimcr 1997; Whitworth et al. 1995). The formation of this 
first star generally occurs in the lowest of the local potential 
minima. Surrounding condensations with slightly lower gas 
densities form additional staxs (e.g. in the filaments whose 
intersection generated the first dense clump). Both the stars 
and the residual gas arc attracted by their mutual gravita- 
tional forces and fall towards each other. The interactions 
between the gas and the protostars dissipate some of the ki- 
netic energy of the latter (Bonnell ct al. 1997), allowing the 
stellar objects to rapidly come close to the initial star and 
its disc-born companions to form a high-density sub-cluster 
containing from 2 up to 8 stars: a small- A'^ cluster. This pro- 
cess repeats itself in other parts of the cloud. Given the size 
of the cloud wc have modelled (and consequently, the num- 
ber of Jeans masses initially present in the system), no more 
than 3 of these star-forming sites are ever produced. Subse- 
quently, sub-clusters are attracted to each other and merge 
to form the final mini-cluster. Thus, the star formation pro- 
cess is hierarchical in nature, as has been vividly illustrated 
(for a 1000 Mq cloud) by Bonnell, Bate & Vine (2003). 

These sub-clusters bear much similarity with the small- 
A'^ clusters modelled by Delgado-Donate, Clarke & Bate 
(2003; henceforth DCB03) in the sense that, initially, the 
cluster components are arranged in a dynamically unstable 
configuration. As the systems seeks to attain stability, the 
cluster breaks up: the most massive components form a hi- 
erarchical multiple (typically a binary or a triple system) 
whilst the low mass components are ejected, either to largo 
separations or from the cloud altogether. In contrast with 
the DCB03 simulations, the number of stellar seeds is not 
fixed and therefore further star formation events can take 
place, introducing dynamical instability again in the system. 
Sub-cluster merging provides an additional complication to 
the simple small- A'^ cluster picture, as a multiple system can 
be driven to the proximity of another one so that further 
interactions take place. 

The fact that low mass components are the prime candi- 
dates to be ejected means that, given the high stellar density 
of each sub-cluster and the tendency of sub-clusters to inter- 
act with each other, few bound pairs involving a low mass 
component - i.e. low mass ratio, wide or low mass pairs 
- can survive the interaction with other cluster members. 
The binding energy of these pairs is too low. In addition, in 
the present simulations accretion discs surround most of the 
stellar objects formed, and when surrounding a binary, tend 
to drive the mass ratio q to high values (Bate & Bonnell 
1997). After sub-cluster merging, a new process comes into 
play: exchange of binary components (Valtonen & Mikkola 
1991). Thus, the lightest companions are exchanged by more 
massive ones and hence, the probability for the surviving 
binaries to have nearly equal-mass components is enhanced 
(see also Bate, Bonnell & Bromm 2002b). In addition, where 
quadruples are formed involving binary-binary pairs, the 
surrounding circum-quadruple disc will tend to drive the 
total masses of each binary to similar values. 

A large fraction of the bodies ejected from the unstable 
sub-clusters become sub-stellar objects. This is so because 
the ejected objects are typically the low mass components of 



the system and because, after being ejected from the dense 
region in which the multiple system sits, the ejected body 
is largely deprived of further accretion (Reipurth & Clarke 
2001; Bate, Bonnell & Bromm 2002a). Hence, taking into 
account that all the objects in these simulations start with 
a mass close to the opacity limit (a few Jupiter masses), 
it is not surprising that many of the ejectae become brown 
dwarfs. 

The spatial distribution of stars and brown dwarfs 
within each multiple (henceforth the configuration) is shown 
in Table 1. For each simulation, the configuration column 
contains two rows: the upper one refers to the internal struc- 
ture of the multiples at t — 0.5 Myr, i.e. at the end of the 
hydrodynamic calculation (hereafter hydrodynarnic configu- 
ration); and the bottom row corresponds tot = 10.5 Myr, i.e 
at the end of the A'^-body integration. All the other columns, 
except the first two, also refer to both the hydrodynamic and 
N-body results. 

3.1.1 Hydrodynamic results 

It can be seen from the hydrodynamic configurations that the 
basic building blocks of the multiples' internal hierarchy are 
binary stars. This property refiects the origin of these bound 
systems as byproducts of a hierarchical formation process, 
in which the typical outcome of a small- A^ cluster disinte- 
gration is a tightly bound binary star. Subsequent merging 
of several small-A'^ clusters bind two or more binaries into 
one single multiple system. Column 3 ([Nm ; Ns]) shows on 
the left the number of stars and brown dwarfs that remain 
in multiple systems at the end of the simulations; in the 
case that more than one mutually unbound multiple system 
is formed, the membership number of each is separated by 
commas within brackets. On the right is shown the number 
of stellar objects that escaped from the cloud. Overall, the 
ratio of unbound to bound objects is approximately 2/3 at 
0.5 Myr. 

At an age of 0.5 Myr, 60% of the stars and brown dwarfs 
are locked in 12 multiple systems, with about a third of the 
components being low-mass, weakly bound outliers. Exclud- 
ing these outliers and unbound singles, 7% of the remaining 
objects are in pure binaries (2 systems), 14% are in quadru- 
ples (2 systems), 35% are in quintuples (4 systems), 32% 
are in sextuples (3 systems) and 12% are in multiples with 
seven components (1 system). 

Figure 1 shows the value of the mass ratio q (left panel) 
and semi-major axis a (right panel) of the resulting multiples 
as a function of primary mass (see figure caption for an ex- 
planation of the symbols) . It can be seen that no binary star 
with primary mass larger than ~ 0.2 Mq is formed, despite 
the initial mass of the components being 200 x smaller. In 
addition, all except one binary have q larger than 0.5. The 
exception appears in the region located between the solid 
and dashed lines, where binaries with brown dwarf compan- 
ions are found. Apart from this binary, brown dwarf com- 
panions are only found as wide components of high-order 
multiples. 

The right panel illustrates two trends. Firstly, high-A'' 
multiples arc typically wider than low-A" multiples. This 
trend reflects the hierarchical or nearly hierarchical spa- 
tial distribution that the multiple's components adopt for 
the system to attain dynamical stability - note that this 
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Figure 1. Mass ratios and semi-major axis (at 0.5 Myr). Left panel: mass ratio q vs primary mass (in Mq), for all the systems formed 
in the simulations. Diamonds correspond to binaries (either independent or bound to larger structures), triangles denote triples (same 
as for binaries), squares quadruples (same), asterisks quintuples (same) and crosses higher-order multiples. The values of q must be 
interpreted as meaning the ratio of the mass of the outermost object of the sub-system under consideration to the total mass of all the 

objects interior to it (i.e. for a triple [(l-2)-3], q = m3/(ml+m2)). Primary mass refers to the total mass of all the objects interior to 
that for which the mass ratio is being calculated. Between the solid and the dashed line are located those systems with a brown dwarf 
companion. Brown dwarf binaries would appear to the left of the dashed line. Right panel: semi-major axis (in AU) vs primary mass (in 
Mq) for all the multiple systems produced by the simulations. Symbols as in left panel. 
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trend is better visualised if the symbols rather than the 
primary masses are considered. The second trend reflects 
a shortcoming of these simulations: binaries have a mean a 
of « 10 — 20 AU, i.e. the vast majority are close binaries 
and only 1 system has a a > 50 AU. This result is inde- 
pendent of primary mass. In other words, wide pure bina- 
ries of any mass axe underproduced by our simulations and 
although closer binaries might have been formed if the soft- 
ening length used to smooth the gravitational force at short 
distances had been chosen smaller, no wider system can be 
produced with these initial conditions. The reason for this 
is twofold: first, there is not enough global angular momen- 
tum initially to form many wide binaries or a very wide 
binary, and second but most importantly, the high density 
of the sub-clusters prevents wide systems from surviving. 
This conclusion can be applied to all the systems underpro- 
duced by our simulations: low mass, low mass ratio and wide 
binaries share the property of having low binding energies. 
Only the systems with large binding energies can survive 
the dynamical encounters taking place in the dense clusters 
formed by collapsing turbulent flows. 

Figure 2 shows pictorial sketches of the spatial distri- 
bution of the components of a representative sample of the 
multiple systems resulting from our simulations, i.e. a bi- 
nary star orbiting a triple (top left), a binary quadruple (a 
binary orbiting another binary) with some outliers (bottom 
left), a binary orbiting a binary quadruple plus some out- 
liers (top right), and a planetary quadruple (bottom right). 
The first three classes of multiples comprise about 75% of all 
the multiples formed in the simulations. The remaining 25% 
(the fourth class) is comprised by planetary multiples, sys- 
tems in which companions are not members of binary/triple 
systems other than the multiple itself. One of the noteworthy 
features of the mass distribution within the multiples is how 



similar the masses of all components are, either taking an 
individual star or a binary star as the fundamental unit. In 
other words, most binaries have very high q and, when two 
binaries are bound to each other, the ratio of the total mass 
of one to the total mass of the other is also close to unity (in 
the range 0.5 — 1; see left panel of Figure 1). This pattern of 
mass distribution can be readily related to the hierarchical 
formation mechanism that produces the multiples: at each 
sub-cluster merging event, the binary stars involved inter- 
act strongly, exchanging components, and accreting material 
with higher specific angular momentum than that of their 
orbit. These two processes invariably favour the formation 
of bound pairs (a pair of stars in the first small-N cluster, a 
pair a binaries after a merging event, etc.) with high values 
of g (> 0.5). 

Whereas the flrst star formation episodes involve as 
much direct fragmentation in filaments Eis disc fragmenta- 
tion (BBB), at later times, star formation is dominated by 
instabilities in circum-binary /quadruple discs. Most objects 
formed at late times are unlikely to remain bound since they 
start off with a mass close to the opacity limit for fragmen- 
tation whereas the other sub-cluster members have been ac- 
creting for some time and have much larger masses, therefore 
being better positioned for future three-body encounters. 
The net result is that only those objects formed during the 
first free-fall times after star formation begins are likely to 
remain in bound structures, whilst the rest will be ejected at 
very large separations {outliers) or simply escape from the 
cloud. 

This is reflected in Table 1 in the penultimate column 
(Nout), which gives the number of objects that being bound 
to a given multiple system, have a mass much lower than 
the other components, and highly eccentric orbits. Outliers 
orbit the CoM of the multiple at large distances, « 10^ — lO"* 
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Figure 2. Pictorial sketches of the spatial distribution of the components of a representative sample of the multiple systems resulting 
from our simulations. Large dots represent stars in the mass range 0.1 — 0.6 Mq and small dots represent outliers (M < 0.1 Mq). 
The number beside each dashed straight line gives the distance (in AU) between the sub-systems joined by that line. These are typical 
configurations at an age of 0.5 Myr. On the bottom right a planetary multiple is shown; the other systems involve two or more mutually 
bound binaries. 




AU (see Figure 1 [right panel; crosses correspond to outliers] 

or the last column of Table 1). The population of outliers 
is mostly comprised by brown dwarfs, and makes up about 
20% of all the objects in bound systems. 

3.1.2 N-body results 

The A''-body evolution of the stellar systems result in the 
break-up of many of the higher-order multiples. At an age 
of 10.5 Myr, 36% of all stars and brown dwarfs are locked in 
18 multiple systems, with only about 2% of the components 
being low-mass outliers. Excluding these, 42% of the remain- 
ing bound objects are in pure binaries (11 systems), 12% are 
in triples (2 systems), 15% are in quadruples (2 systems), 
19% arc in quintuples (2 systems) and 12% are in sextuplcs 
(1 system). The companion frequency /c or average number 
of companions per star system is defined as ' 
where s represents the number of singles, b the number of 
binaries, t the number of triples, q the number of quadru- 
ples and k the number of quintuples. This quantity changes 
substantially during the 10 Myr of A/'-body evolution, from 
1.0 to 0.3. On the other hand, the multiplicity frequency /m, 
defined as , shows just a small variation, from 



0.18 to 0.17. The multiplicity frequency shows such small 

variation because the number of multiples and the number 
of singles increase during the A''-body evolution by similar 
factors, « 1.4 and k, 1.5 respectively. Since these factors 
appear one in the numerator and the other in the denom- 
inator of /m, they cancel out almost completely and thus 
/m remains practically unchanged. These two quantities, /m 
and predominantly /c, refiect the process of multiple disinte- 
gration that has taken place during the A-body evolution: a 
substantial number of high-order multiples have broken-up 
into their fundamental units, i.e. pure binaries, which now 
comprise the largest group. The fraction of multiples out 
of the total number of systems remains almost constant, 
whereas the average number of members per system de- 
creases by a factor of 3. This decrement is predominantly 
due to the ejection of outliers (only 3 survive after 10 Myr), 
in other words, the transference of bound objects to the pop- 
ulation of singles. This process is responsible for « 80% of 
the variation in /c. Ejections of outliers have a side eflfect: 
they harden the inner multiple, which very often is a binary 
quadruple. After several ejections, the separation between 
the binaries may have diminished enough so as to render 
the multiple unstable. After several orbital timescales, close 
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binary-binary interactions result in the disintegration of the 
multiple, cither into two separate binaries or a binary and 
two singles. This process, however, only accounts for « 20% 
of the variation in /c, which is clearly dominated by the 
outlior-to-singlc transference. A noteworthy feature of this 
disintegration process is that the bottom level of the multi- 
ple's hierarchical distribution, namely the binary stars, re- 
mains substantially intact. For example, in simulation a3C 
two binaries and one triple are mutually bound at 0.5 Myr. 
After Af-body integration to 10.5 Myr, these systems have 
become mutually unbound, but their internal configuration 
is still the same. There is only one case {a3D) in which the 
disintegration of the multiple results in the exchange of com- 
ponents between two binaries together with the disruption 
of one of them. 

The typical orbital timescale of the outliers is ~ 10® yr. 
After a few orbits these objects should experience several 
close encounters with the inner multiple (since they have 
very eccentric orbits), thus making ejection likely to happen 
during the next few x 10® yr after the start of the A^-body 
integration. It is not surprising then that at an age of 10.5 
Myr only three out of the original 30 outliers have remained 
in bound orbits, and a substantial fraction of the multiples 
have disintegrated. Overall, the ratio of unbound to bound 
objects is approximately 5/3 at 10.5 Myr. In addition, the 
number of single brown dwarfs has increased by w 60% in 
the same period. 

We find that the ratio of binary to higher-order mul- 
tiples at 0.5 Myr is 2/11, but it rises to 10/8 at 10.5 Myr. 
DM91 find a ratio binary:higher-order multiple of « 5 : 1 
for the solar neighbourhood. However, their sample might 
suffer from a bias against N > 2 systems; e.g. some of the 
systems might be triples disguised as binaries either because 
one of the components might be itself an unresolved binary 
or because a faint tertiary companion might remain unde- 
tected. DM91 did spectroscopic follow-up of stars but the 
wider systems were drawn from the literature; there could 
be many missing wide components, which in turn might 
themselves be multiples. In fact. Mayor & Mazeh (1987) 
and Mazeh (1990) suggest that 25-50% of the spectroscopic 
binaries might be higher-order multiples. The evidence for 
this is deduced from the study of some tight binaries which 
exhibit non-negligible residuals after their orbit is fitted by 
a Keplerian 2-body orbit. Moreover, spectroscopic studies 
by Tokovinin and coworkers (Tokovinin 1997; Tokovinin & 
Smekhov 2002) suggest that a substantial fraction of the 
components of binary or triple systems might be themselves 
unresolved binaries. Therefore, it seems unclear at the mo- 
ment whether our results are in agreement with the observed 
ratio of binary:higher-order multiples, as deduced from the 
studies on individual systems. However, as we will show 
later, our results compare well with the photometric width 
of the main sequence of young clusters like Praesepe. There- 
fore, the high incidence of A > 2 multiples that we predict, 
given the present observational constraints, cannot be ruled 
out. 

At the end of the AT-body integration 17 out of the 
18 multiples have stable configurations (for at least several 
xlO* yr), according to the criteria of Egglcton & Kisclcva 
(1995), whilst the 18*'' may further decay on a timescale 
of a few Myr. However, external perturbations may also 
affect the stability of these systems. It is expected that a 



substantial fraction of all multiple systems ever formed out 
of a cloud core will become members of a larger structure, 
namely a young stellar cluster or association, or perhaps an 
open cluster in the long term (Clarke et al. 2000, Lada & 
Lada 2003). Dynamical interactions within these larger stel- 
lar aggregates may produce the disruption of some of the 
widest systems that survived the AT-body integration (e.g. 
as in the Orion Nebula Cluster; Scally, Clarke & McCaugh- 
rean 1999). In a longer term still, for those systems joining 
the field population, the gravitational field of the galaxy will 
act to disrupt wide systems (Weinberg, Shapiro & Wasser- 
man 1987), limiting the maximum separation between com- 
ponents to a value somewhere around 10"* AU for solar-mass 
stars (Close et al. 2003). 

Single and binary stars attain comparable velocities in 
the range 1 — 10km a~^. Higher-order multiples display lower 
velocity dispersions. This kinematic segregation as a func- 
tion of A is the expected outcome of the break-up of unsta- 
ble multiples, whereby the ejected objects (typically singles, 
or less often binaries) acquire large velocities whereas the 
remaining more massive multiple recoils with a lower speed. 
Therefore, we would expect low mass SFRs like Taurus or 
Ophiuclms to display an overabundance of A^ > 2 systems 
as lower- Af systems may escape more easily the potential 
well of these associations. A more detailed discussion on the 
kinematic properties of the stellar objects resulting from our 
simulations is given by DCB04. 

3.2 Multiplicity fractions 

Overall, the multiplicity fraction /m derived from these sim- 
ulations is low: it takes a value close to 0.2, much lower than 
that observed in clusters and the field. However, the variar 
tion of the multiplicity fraction with primary mass can be 
very steep, as Sterzik & Durisen (2003) have shown, both 
observationally and theoretically, and a low total value may 
thus mask high and low values in different mass ranges. Fig- 
ure 3 shows the dependence of /m on primary mass, for 
the hydrodynamic simulations exclusively, at 60% efficiency 
(squares joined by a solid line) and 30% efficiency (circles 
joined by a dashed line). We pointed out in Section 3.1 that 
the overall value of /m at 0.5 Myr is very similar to that 
found at 10.5 Myr. This conclusion also applies to the indi- 
vidual values of /m in each mass bin. Therefore, the 0.5 Myr 
60% efficiency curve is also a good representation of the /m 
dependence on primary mass 10 Myr later. Also shown in 
Figure 4 are the observational results from DM91, Fischer 
& Marcy (1992; FM92), Marchal et al. (2002, M&:02) and 
Bouy et al. (2003), Close et al. (2003), Gizis et al. (2003) and 
Martin et al. (2003) [BCGM03]. As expected, the multiplic- 
ity fraction is an increasing function of primary mass, with 
the most massive stars produced in these simulations hav- 
ing a multiplicity fraction very close to 1. At the other end, 
brown dwarfs are rarely binary primaries. Nevertheless, the 
shape of the /m curve is sensitive to the efficiency assumed. 
For 30% efficiency, the results show better agreement with 
observations at the sub-stellar regime, while the 60% effi- 
ciency results match more closely the observed multiplicity 
fractions for KM stars. 

Two points must be stressed however: first, the limita- 
tions imposed by the initial total mass in gas of the cloud, set 
to 5 Mq. The fact that about 10 stars/brown dwarfs form. 
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Figure 3. Multiplicity fraction /m as a function of primary mass 
(in Mq). Dotted vertical lines separate the different mass bins 
for which fm has been calculated. Squares joined by a solid line 
and circles joined by a dashed line denote the values of fm at 60% 

and 30% efficiency, respectively. The results shown in this plot 
correspond to an age of 0.5 Myr (see text). 
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on average, in each calculation, and that those with higher 
mass will be invariably in bound systems, place a /m value 
close to 1 in the solar-mass stars bin. For low-mass star- 
forming regions (SFR) such as Taurus (Duchene et al. 1999) 
or p Ophiuchus (Barsony, Koresko & Matthews 2003), our 
choice of 5 Mq clouds seem to be appropriate: both SFRs 
lack stars more massive than about 2 M0 and similarly have 
a binary fraction approximately twice as high as that of field 
stars (which stands at « 50%), very close to what we find. 
For more massive SFR, models of more massive clouds are 
more appropriate: had we modelled a cloud with 50 or 500 
M© , the highest value of the multiplicity fraction would have 
moved to higher primary masses (e.g. 5 to 10 M©; Bonnell 
& Bate 2002). This would bring closer agreement with the 
results in the field (DM91; Shatsky & Tokovinin 2002). 

Second, the discrepancy with observations at the low 
mass end (too few brown dwarf binaries) cannot be allevi- 
ated by modelling more massive clouds. BBB, modelling a 50 
M0 cloud, found only one brown dwarf binary candidate out 
of ~ 20 single brown dwarfs, giving them a 5% binary frac- 
tion at the sub-stellar regime. In practice, this binary would 
have likely been disrupted or simply accreted to stellar val- 
ues had the simulation been carried forward longer, since we 
do form many of these binaries in our simulations but they 
soon are destroyed or become stellar binaries. It remains a 
puzzle how to produce a fraction of very low-mass or brown 
dwarf binaries as high as is observed. Certainly, clouds with 



an initial mass of ~ 0.5 M0 and a lower number of Jeans 
masses than used in the present models, do produce many 
of these very low-mass binaries. However, in real clouds, one 
would expect very low mass cores to interact strongly with 
others, as in these simulations, and only the prompt ejec- 
tion of a brown dwarf binary via an encounter with a more 
massive binary could result in the formation of these sys- 
tems. However, neither our models nor those of Bate, Bon- 
nell & Bromm (2002a) seem to favour this scenario: ejection 
of low mass binaries, as opposed to low mass singles, are 
rare. One possibility is that the resolution limit imposed by 
our choice of the softening radius of the sink particles may 
overestimate the number of disruptions taking place among 
low-mass binaries. Binaries with separations smaller than 
« 1 AU cannot form, since gravitational accelerations are 
smoothed within that radius. Thus, there is a limit to the 
minimum binding energy that a bound pair can have. The 
lower the binding energy of a low-mass bound pair, the like- 
lier it is its survival to encounters with higher mass stars. If 
a brown dwarf binary could become tight enough (i.e. have a 
separation < 1 AU), it might be able to survive as a bound 
pair an encounter and be ejected. This possibility will be 
explored in future simulations. 



4 OBSERVATIONAL TESTS 

Some of the properties of our simulated multiple systems can 
be immediately compared with current observational data, 
thus providing an useful benchmark against which to assess 
the validity of our models. Some other properties, however, 
do not seem to have been looked upon extensively by ob- 
servers and may require additional observational measure- 
ments. Therefore, in addition to the conventional cross-check 
with observations, the aim of the next subsections is also to 
sketch a series of observational tests that, if undertaken, 
might help to discriminate among different theoretical mod- 
els of low-mass star and brown dwarf formation. 

4.1 Colour-magnitude diagrams 

Colour-magnitude (CM) diagrams provide useful informa- 
tion on the binarity of a given homogeneous sample of stars 
and brown dwarfs. Two different main sequences are fre- 
quently discovered in these diagrams, one corresponding to 
single stars or individual stars in wide (i.e. larger than the 
instrument resolution limit) binary systems, and the other 
usually ascribed to unresolved binary stars. If more than 
two stars contribute significantly to the total luminosity of 
a multiple system, the corresponding point on the CM di- 
agram might be expected to have a lower magnitude than 
that of a pure binary of the same colour, thus broadening the 
binary sequence. Our simulations predict the existence of a 
substantial number (50% at 10.5 Myr) of multiple systems 
made up of three or more stars with comparable masses. 
Thus, we would like to know whether a CM diagram of our 
simulated stellar systems results in a binary sequence whose 
width is consistent with those of real systems or not. 

Figure 4 shows the CM diagram corresponding to our 
multiple systems at 10.5 Myr, after converting masses to / 
and K magnitudes using the tracks by Baraffe et al. (1998). 
We have assumed an age of 600 Myr for the tracks and a 
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resolution of 200 AU (components of a multiple system sepa- 
rated by more than 200 AU are considered independent sys- 
tems when computing their magnitudes). All systems closer 
than 200 AU are stable, thus we do not expect any further 
evolution of this diagram for an older system. Red squares 
correspond to unresolved multiple systems, and blue open 
circles to individual stars. Over-plotted, in solid black cir- 
cles, are the infrared observational results from the Prae- 
sepe cluster (« 600 Myr old and located at ~ 180 pc), by 
Hodgkin et al. (1999). The stars used in this diagram have 
been selected from the Hambly et al. (1995) proper-motion 
and photometric survey of Praesepe. The magnitudes of this 
sample have been measured from photographic plates, and 
the resolution limit ranges between 2 "and 3 ", hence our 
criterion of 2.5 "{^ 200 AU at 180 pc) used to identify a 
given simulated multiple either as unresolved or as a set of 
independent objects. 

Three features from Figure 4 may be highlighted: first, 
a binary sequence is apparent in the simulated data, and its 
width is not significantly different to that of the Praesepe 
cluster, except for systems redder than I — K = 2.5. This 
seems to suggest that the formation of triple, quadruple and 
higher-order multiples of the sort produced by our models is 
not ruled out by observations and indeed might be common 
in real clusters. 

Second, our models fail to produce as many very low- 
mass binaries as are observed. The observed binary sequence 
for Praesepe (as e.g. for the Pleiades; Pinfield et al. 2003) 
does continue well into the sub-stellar regime. As explained 
in Section 3.2, this observational result is unaccounted for by 
our simulations at 60% star-formation efficiency, and repro- 
duced only to some extent if very low efficiencies (probably 
unphysical for a molecular cloud core) are assumed. 

Third, the sequence of simulated singles does not extend 
to colours as blue as those of the multiples. The reason for 
this is that most of the mass of the cloud ends up in the inner 
components of multiple systems (singles are mostly ejected 
objects, which remain bound at large separations or escape 
completely from the cloud and are thereafter deprived from 
further accretion). Had we modelled a more massive cloud, 
ejections would have also occurred among higher-mass ob- 
jects and the singles sequence might have extended to bluer 
colours. Likewise, the lack of wide pure binaries in our results 
contributes to the paucity of objects in the upper part of the 
singles sequence. This latter problem might be solved by in- 
creasing the net angular momentum in our models, which 
initially corresponds to a /3 parameter (the ratio of the rota- 
tional energy to the magnitude of the gravitational energy 
of the cloud) of ~ 10"^ 

4.2 High-resolution imaging 

High-resolution imaging with HST or ground-based tele- 
scopes using adaptive optics has proved very successful in 
resolving tight binary systems. Most recent studies on bi- 
narity among previously unresolved systems have focused 
on low-mass stars and brown dwarfs. Bouy et al. (2003), 
Close et al. (2003) and Gizis et al. (2003) provide the most 
complete set of observations to date of field very low-mass 
and brown dwarf binaries. Alternatively, Martin et al. (2003) 
have surveyed young stellar clusters, such as a Persei and 
the Pleiades, to pin down the binary fraction among the 



Figure 4. Colour-magnitude (CM) diagram (I vs I — K) for the 
stars and brown dwarfs formed in the simulations. Red squares 
represent unresolved multiples and blue open circles singles. Ob- 
servational measurements from the Praesepe cluster are shown as 
black filled circles. Resolution at 200 AU (or 2.5 "at 180 pc) 
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brown dwarf population. All these studies suggest that the 
binary fraction among very low-mass stars and brown dwarfs 
is much lower than for G-type stars, a trend that we also 
find in our simulations (see Figure 3). Specifically, Close et 
al. (2003) find a binary fraction in the range 15% ± 7% for 
M8.0-L0.5 stars with separations greater than 2.6 AU, and a 
semi-major axis distribution peak at ~ 4 AU. In their sam- 
ple, no very- low mass binary has a separation greater than 
15 AU. These results are in agreement with those of Gizis et 
al. (2003) and Bouy et al. (2003). Martin et al (2003) also 
concludes that a binary frequency in the range 10% — 15%, 
a bias to separations smaller than about 15 AU, and a ten- 
dency towards high mass ratios (g ^ 0.7) are the fundamen- 
tal properties of brown dwarf binaries. It must be noted, 
however, that the 10% — 15% observed binary frequency is 
most likely a lower limit and that, in general, there is con- 
siderable uncertainty about the properties of brown dwarf 
binaries. In our simulations, most pure binary systems are 
tight (see Figure 1, right panel), with a median semi-major 
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axis of 10 AU, but few of these systems have primary masses 
below 0.3 M0 (for a discussion of this issue, see section 3.2). 

Adaptive optics imaging surveys of common associa- 
tions, such as Tucana-Horologiurn (Chauvin et al. 2003), 
TW Hydrae (Zuckerman et al. 2001b), MBM 12 (Hearty et 
al. 2000) or the (3 Pictoris moving group (Zuckerman et al. 
2001a) are ideal to study the multiplicity properties of young 
stars, due to their proximity to Earth (< 70 pc) and their in- 
ferred youth (e.g. ~ 20 Myr for the /3 Pictoris group). These 
associations clearly demonstrate the formation of stars in a 
clustered mode and, therefore, represent an invaluable lab- 
oratory against which to compare the results obtained by 
models of clustered star formation. If dynamical interactions 
among embedded stars arc indeed as important as shown 
by our models, this should be reflected most clearly on the 
multiplicity properties of the stars populating these groups: 
i.e. the binary fraction should be higher than in the field 
or in clusters (since many singles should have escaped the 
group), low-mass objects should be found to be effectively 
bound to more massive binaries/multiples, but at large sepa- 
rations, and conversely, single unbound brown dwarfs should 
be rare. Companion frequencies should also be high (in the 
range 0.5 — 1). Recent results by Brandeker, Jayawardhana 
& Najita (2003) point to this direction: they find that the 
multiplicity frequencies in the TW Hydrae and MBM 12 
groups are high (0.58 ± 0.12 and 0.64 ± 0.16, respectively) 
in comparison with other young regions such as Taurus or p 
Ophiuchus (Duchone 1999). The companion frequencies (or 
average number of companions per star system) for these two 
moving groups are even higher: 0.84 ± 0.22 and 0.91 ± 0.30, 
suggesting that many of the multiples have N > 2. These 
results may be compared with the value of /c that our mod- 
els predict if e.g. 50% of the singles have left the group at 
t = 10.5 Myr, i.e., /c = 0.7. 

High-resolution imaging has recently provided the first 
two brown dwarf candidates orbiting solar-type stars at 
small separations (Els et al. 2001; Liu et al. 2002). These 
brown dwarfs are at distance of « 15-20 AU of the primary. 
Although these observations imply that brown dwarf com- 
panions do exist at separations comparable to those of giant 
planets in our solar system, the frequency of this type of sys- 
tems remains unclear: half of the stars harbouring planets 
have long term trends in their radial velocities due to unseen 
companions (Fischer et al. 2001), but so far no other brown 
dwarf companion at close separations heis been found. In our 
simulations we find one binary (Table 1, ref. (6 — 7) in a3C) 
out of 24 (at 10.5 Myr) to have a close brown dwarf compan- 
ion: this gives a frequency of ~ 4%. This binary is composed 
of a 0.21 Mq star and a 0.05 M© brown dwarf {q = 0.24), 
has a semi-major axis of 10 AU and an eccentricity of 0.7. 

4.3 Spectroscopy 

Radial velocity surveys find that the incidence of brown 
dwarfs within 4 AU of solar- type FGK stars is < 1%. The 
evidence for this brown dwarf desert at very small separa- 
tions first emerged from radial velocity surveys in the late 
80's and early 90's (e.g Walker et al. 1995), and has been 
confirmed by current high-precision radial velocity programs 
(Marcy & Butler 2000; Halbwachs et al. 2000). Our results 
are consistent with this observed brown dwarf desert: no bi- 
nary has a sub-stellar companion closer than 10 AU. The 



explanation of this result is very simple: most brown dwarfs 
companions to stars are formed via the fragmentation of cir- 
cumstellar discs, but this same disc drives quickly the mass 
ratio of the binary towards unity, due to the accretion of 
material with higher specific angular momentum than that 
of the binary orbit (Bate & Bonnell 1997; Bate 2000). In 
some other cases, the brown dwarf companion is also likely 
to be exchanged by a more massive object after a three-body 
encounter. 

Radial velocity measurements of the components of vi- 
sual double and multiple stars (Fekel 1981; Tokovinin 1997; 
Mazeh et al. 2001; Tokovinin & Smekhov 2002) have resulted 
in the detection of a substantial fraction of spectroscopic 
sub-systems. Tokovinin (1997) points out that the frequency 
of sub-systems among spectroscopic binaries with periods 
shorter than 10 days is at least of 40%. More recently, 
Tokovinin & Smekhov (2002) have found a frequency of spec- 
troscopic sub-systems of « 20% per component among 26 
resolved visual binaries (i.e. 20% of these binaries turned out 
to be triples or alternatively, 10% may be binary quadru- 
ples). In addition, 18 out of 59 apparently single tertiary 
components turned out to be spectroscopic binaries (i.e. 30% 
of the systems are quadruples instead of triples). These re- 
sults seem to be in agreement with the high abundance of 
high-order multiple systems found in our simulations. Note 
that although no binaries with semi-major axis smaller than 
1 AU can be produced by our simulations (due to soften- 
ing of the gravitational acceleration between point masses 
at short distances) there is no reason to believe that such 
spectroscopic binaries would not have formed had the reso- 
lution limit been much lower; therefore, a direct comparison 
of Tokovinin & Smekhov findings with our results is indeed 
meaningful. Further radial velocity observations of the com- 
ponents of multiple systems would be very helpful to pin 
down the real fractions of binary, triple and quadruple sys- 
tems, and thus constrain different models of star formation. 

4.4 Very low-mass stars and brown dwarfs at 
large separations 

In Section 3.1 we pointed out that, in our 0.5 Myr mod- 
els, very low mass objects (often brown dwarfs) are very 
common as companions to more massive stars, but at large 
separations. Furthermore, most of these outliers are bound 
not to single stars but to close binaries, triples or binary 
quadruples. At an age of 10.5 Myr, the fraction of multi- 
ple systems exhibiting brown dwarfs at large separations 
is much lower: 3 out of 18, in contrast with the result of 
10 out of 13 found at 0.5 Myr. Nevertheless, 3 out of the 
4 brown dwarfs that remain bound do so orbiting binaries 
or higher-order multiples at large separations. The excep- 
tion is the brown dwarf secondary found in simulation a3C, 
which is only 10 AU away from the primary. Therefore, it 
appears that many bound brown dwarfs (~ 3/4), and most 
bound brown dwarfs in wide orbits (~ 3/3), should be orbit- 
ing binary, triple or quadruple systems. Currently, 12 brown 
dwarf/very low mass companions at wide separations are 
known (Kirkpatrick et al. 2001a,b; Wilson et al. 2001), prov- 
ing that the brown dwarf desert docs not extend to large 
separations (Gizis et al. 2001). Two of these brown dwarfs 
(G1337C and G1584C) are orbiting visual binaries, at w 900 
and 3600 AU respectively, and a third (G1570D) is bound to 
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a triple system (at a distance of 1500 AU). That is, 3 out of 
the 11 brown dwarfs orbiting stars at large separations are 
known to orbit a binary or triple system. A clear prediction 
of our simulations can be applied to the other 8 systems in 
which a brown dwarf is apparently orbiting a single star. A 
large fraction of these singles should turn out, after closer 
examination, to be iV > 2 multiples (i.e. a spectroscopic bi- 
nary, triple, etc.). This could be tested very simply, since a 
significant fraction of all spectroscopic binaries are known 
to be twins - i.e. have nearly equal-mass components (Halb- 
wachs et al. 2003) -, as is also the case in our models, and 
twins are easier to detect spectroscopically than low mass 
ratio binaries. An observational case in line with our pre- 
dictions has been described by Brandckcr, Jayawardhana 
& Najita (2003) who have recently shown that the brown 
dwarf TWA 5 B is bound at « 120 AU to TWA 5 A, which 
in turn is resolved into a very tight, 3 AU separation, binary 
(or possibly a triple; Mohanty, Jayawardhana & Barrado y 
Navascues 2003). 



5 CONCLUSIONS 

We have undertaken the first hydrodynamical -|- A/^-body 
simulations of multiple star formation that have produced 
a statistically significant number of stable hierarchical sys- 
tems, with component separations in the range ~ 1 — 1000 
AU. These simulations have demonstrated that multiple star 
formation is a major channel for star formation in turbu- 
lent flows. The hydrodynamical simulations are followed for 
~ 0.5 Myr; subsequently, the remaining gas is removed and 
the stellar systems followed as A^-body ensembles for an ad- 
ditional 10 Myr. At this point, all but one of the surviv- 
ing multiple systems are stable, according to the criteria of 
Eggleton & Kiseleva (1995). We find that the properties of 
the resulting multiple systems arc not significantly sensitive 
to the large scale geometry of the cloud - determined by 
the turbulence - but rather to the dynamical and competi- 
tive accretion processes taking place within the mini-clusters 
formed out of the collapse and fragmentation of the cloud. 

At an age of 0.5 Myr, we find that about 60% of stars 
and brown dwarfs are in multiple systems, with about a third 
of these being low mass, weakly bound outliers. Excluding 
these outliers and unbound objects, 7% of the remaining ob- 
jects are in pure binaries (2 systems), 14% are in quadruples 
(2 systems), 35% are in quintuples (4 systems), 32% are in 
sextuples (3 systems) and 12% are in multiples with seven 
components (1 system). The companion frequency is there- 
fore very high, « 1. We find that our multiples consist of 
hierarchies of binaries and triples and that planetary multi- 
ples (in which companions are not members of binary /triple 
systems other than the multiple itself) are comparatively 
rare (occurring ~ 25% of the time). There is a distinctive 
pattern of mass distribution within these multiples, with the 
mass ratio within binaries, and the mass ratios between bi- 
naries, rarely deviating far from unity (values of 0.5 — 1 are 
typical). On the other hand, such systems are typically or- 
bited by several low mass outliers (typically at separations 
of - 10* AU) on eccentric orbits. About 90% of these ob- 
jects are unstable in tirnescales of a few x 10® yr (i.e. a few 
X their typical orbital timescale). 

We find that the 40% of objects that are unbound are 



overwhelmingly of low mass (median mass ~ 0.02 M0). Thus 
our results imply that in the stellar regime, most stars are in 
multiples (« 80%) and that this multiplicity fraction /m is 
an increasing function of mass. In this latter respect, these 
results are qualitatively consistent with a large body of pre- 
vious works on the decay of small- A?" systems, both with and 
without gas (van Albada 1968; McDonald & Clarke 1993, 
1995; Sterzik & Durisen 1998, 2003; DCB03). The high 
values for GK stars are consistent with adaptive optics mea- 
surements of nearby young associations such as MBM 12 
and TW Hydrae (e.g. Brandeker, Jayawardhana & Najita 
2003), where multiplicity fractions as high 0.64 are found, 
and radial velocity surveys of visual binaries (e.g. Tokovinin 
& Smekhov 2002) which raise the percentage of spectro- 
scopic sub-systems to at least 40%. Low- mass SFR such as 
Taurus or p Ophiuchus also show companion frequencies in 
the range 0.3 — 0.5, comparable to those predicted by our 
models at later times. It must be pointed out that the val- 
ues of the multiplicity fraction /m for each mass range do 
not change significantly during the TV-body evolution of the 
systems. 

At an age of 10.5 Myr the fraction of bound and un- 
bound objects has reversed: 40% remain in multiples and 
60% are singles. The companion frequency has dropped to 
Ri 0.3 due to the ejection of bound outliers to the field. This 
transference of objects from bound to unbound orbits results 
in an increase of the number of free floating brown dwarfs 
by « 60%. In this 10 Myr time-span, many multiple sys- 
tems also experience internal decay: excluding the remain- 
ing 3 outliers, 42% of the remaining bound objects are in 
pure binaries (11 systems), 12% are in triples (2 systems), 
15% are in quadruples (2 systems), 19% are in quintuples (2 
systems) and 12% are in sextuples (1 system). 

We pointed out that low mass stars (and, especially, 
brown dwarfs) are locked up in multiple systems at early 
times and subsequently released into the field. This remark 
needs some qualification however. Multiple star formation is 
hierarchical in our simulations, with structures forming on 
a particular scale being modified as a result of subsequent 
merging with structures on a larger scale. It is notable that 
we find low mass outliers at a separation that is similar 
to the initial size of the core, and we speculate that if we 
had modelled a larger volume of cloud, rather than isolated 
cores, we might have found that these outliers would already 
have been stripped off by interactions with structures on 
a larger scale before some of them could settle into stable 
orbits. This suggests that such outliers should be sought 
in imaging of relatively isolated pre-main sequence groups 
and may explain why no such outliers have been detected 
through deep imaging of multiple systems in Taurus (G. 
Duchene, private communication). 

If such outliers do indeed survive the formation process, 
then about 10% of them are in stable hierarchical orbits at 
10.5 Myr. We would thus expect some brown dwarfs to re- 
main in the outer reaches of multiple systems even in the 
field. Our simulations can thus accommodate the existence 
of systems with brown dwarfs as wide companions (Gizis et 
al. 2001) but only if the primaries of these systems are them- 
selves multiple systems. (We find that three out of the four 
bound brown dwarfs present at t = 10.5 Myr are orbiting 
a multiple system). We therefore predict that the primaries 
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of binaries containing a brown dwarf in wide orbit should 
themselves be multiple systems. 

We have examined how well the products of our simula- 
tions compare with the properties of real stellar systems as 
deduced from the colour magnitude diagram of young clus- 
ters (specifically the infrared colour magnitude diagram for 
Praesepe). Our simulation data compares very favourably 
with the width of the main sequence in the mass range 0.4-1 
Mq; indeed the spread of the main sequence in this mass 
range appears to require that stars are commonly assem- 
bled into high order multiple systems, although the number 
of outliers from a pure equal mass binary sequence is not 
large. 

The comparison with the Praesepe colour magnitude 
diagram however illustrates two problems with the simula- 
tion results. Firstly, the simulation produces no single star 
sequence at masses greater than 0.35 M© (colour bluer than 
I — K — 2; there are only two blue circles above 0.35 Mq, at 
0.581 and 0.584 M0), whereas the observational data shows 
such a sequence, indicating that single stars and/or low mass 
ratio binaries are produced in this colour range. We could 
probably alleviate this problem by modelling a more mas- 
sive cloud. This would increase the maximum mass of stars 
produced and enable some stars in the colour range con- 
sidered to be ejected as singles by encounters with more 
massive stars. Although this would bring closer agreement 
with the observed colour magnitude diagram, the lack of low 
mass ratio binaries in our simulations is in conflict with in- 
dependent evidence from field binary surveys, such as that 
of DM91 for G dwarfs. The DM91 data (containing bina- 
ries with median separation 30 AU) showed a distribution 
with mass ratio (g) that rose with decreasing q towards their 
completeness limit of g ~ 0.2, whereas our mass ratio dis- 
tribution (with the exception of the very low mass outliers) 
is strongly concentrated between 0.5 and unity. This inabil- 
ity to reproduce enough extreme mass ratio systems is a 
feature of all hydrodynamic modelling of multiple systems 
to date. It has not however been widely discussed previ- 
ously, since BBB focused on the close systems, where the 
observed mass ratio distribution is in any case much more 
concentrated towards unit mass ratio than for the binary 
population as a whole (Mazeh et al. 1992; Halbwachs et al. 
2003) and where the simulations very naturally reproduce 
the observed brown dwarf desert at very small separations 
(Marcy & Butler 2000; Halbwachs et al. 2000). 

The second problem revealed by the colour magnitude 
diagram relates to the fact that the simulations produce 
almost no binaries with primaries redder than / — 7^ = 2.5 
whereas the data exhibits a pronounced scatter (consistent 
with large numbers of binaries) in this mass range (0.15 — 
0.4Mq). This scarcity of binaries with low mass primaries in 
our simulations also conflicts somewhat with the results of 
binary surveys among M dwarfs (Fischer & Marcy 1992; see 
Figure 4, this paper) and brown dwarfs (Close et al. 2003; 
Gizis et al. 2003; Martin et al. 2003; Bony et al. 2003). 

In summary, we have found that our simulations pro- 
duce large numbers of hierarchical multiple systems and 
that relatively isolated young multiples may harbour weakly 
bound brown dwarf outhers, as a relic of the hierarchical 
formation process in turbulent flows. We predict that where 
brown dwarfs are found in wide orbits, the primary should 
itself turn out to be a multiple. The simulations are con- 



sistent with a number of observational constraints: the high 
(but presently poorly constrained) incidence of hierarchical 
multiples among field stars and pre-main sequence stars, the 
absence of brown dwarfs as close companions to normal stars 
(the brown dwarf desert) and, at a qualitative level at any 
rate, the positive dependence of the binary fraction on pri- 
mary mass. There are two areas in which the simulations are 
not able to replicate observations properly: the simulations 
under-produce the binary fraction at low masses (M dwarfs 
and brown dwarfs) and also do not generate enough wide 
stellar pairs with low mass ratios. ? 
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